Fungal-clocks: Ophio v. Beau

# SAMPLE NAME
## specify sample name
sample.name <- c("beau", "ophio_cflo", "ophio_kim")
## color scheme for the samples
col.scheme <- c("#5A829F", "#AD212F", "black")

# SCRIPT NAME
## specify the name of the script (folder) where figures will be saved
script.name <- "01_comparing_gene_exp_ophio_beau"

# eJTK OUTPUT
## Set GammaP threshold below which genes are classified as rhythmic
gamma.pval = 0.05
## Set false discovery rate for functional enrichment analyses
FDR = 5

# LOAD DATABASES (TC7)
# 1. TC6_ejtk.db
# Desc: This database contains all ejtk-output for TC6
ejtk.db <- dbConnect(RSQLite::SQLite(), paste0(path_to_repo,"/data/databases/TC6_fungal_ejtk.db"))
# which tables are in the database
src_dbi(ejtk.db)

## src:  sqlite 3.29.0 [/Users/biplabendudas/Documents/GitHub/Das_et_al_2022a/data/databases/TC6_fungal_ejtk.db]
## tbls: beau_rhythmic_genes_12h, beau_rhythmic_genes_24h, beau_zscores_08h,
##   beau_zscores_12h, beau_zscores_24h, ophio_cflo_rhythmic_genes_12h,
##   ophio_cflo_rhythmic_genes_24h, ophio_cflo_zscores_08h,
##   ophio_cflo_zscores_12h, ophio_cflo_zscores_24h, ophio_kim_DD_zscores_24h,
##   ophio_kim_LD_rhythmic_genes_24h, ophio_kim_LD_zscores_24h

#
# 2. TC6_data.db
data.db <- dbConnect(RSQLite::SQLite(), paste0(path_to_repo,"/data/databases/TC6_fungal_data.db"))
src_dbi(data.db)

## src:  sqlite 3.29.0 [/Users/biplabendudas/Documents/GitHub/Das_et_al_2022a/data/databases/TC6_fungal_data.db]
## tbls: beau_expressed_genes, beau_fpkm, beau_log2fpkm, beau_zscores,
##   ophio_cflo_expressed_genes, ophio_cflo_fpkm, ophio_cflo_log2fpkm,
##   ophio_cflo_zscores, ophio_kim_DD_expressed_genes, ophio_kim_DD_fpkm,
##   ophio_kim_DD_log2fpkm, ophio_kim_DD_zscores, ophio_kim_expressed_genes,
##   ophio_kim_fpkm, ophio_kim_log2fpkm, ophio_kim_zscores

#

Overview/Goals

1. General patterns of gene expression

# number of all genes
all.genes <- list()
for (i in 1:length(sample.name)) {
  all.genes[[i]] <- tbl(data.db, paste0(sample.name[[i]] ,"_fpkm")) %>%  
    collect()
  
  writeLines(paste("Number of genes in", sample.name[[i]], ":", nrow(all.genes[[i]])))
}

## Number of genes in beau : 10364
## Number of genes in ophio_cflo : 7455
## Number of genes in ophio_kim : 8577

# A1: genes that have NO expression (FPKM == 0 at all time points)
not.expressed <- list()
for (i in 1:length(sample.name)) {
  not.expressed[[i]] <-
    tbl(data.db, paste0(sample.name[[i]] ,"_fpkm")) %>% 
    collect() %>% 
    filter_at(vars(starts_with("Z")), all_vars(. == 0)) %>%
    pull(gene_name)
  
  # How many genes are not expressed?
  writeLines(paste("n(genes-NOT-EXPRESSED) in", sample.name[[i]], ":", length(not.expressed[[i]])))
  
}

## n(genes-NOT-EXPRESSED) in beau : 759
## n(genes-NOT-EXPRESSED) in ophio_cflo : 190
## n(genes-NOT-EXPRESSED) in ophio_kim : 111

# A2: run enrichment (make plot of enrichment found of non-expressed genes)
for (i in 1:length(sample.name)) {
  writeLines(paste("running GO enrichment for NOT-EXPRESSED genes in", sample.name[[i]]))
  # run enrichment
  not.expressed[[i]] %>% 
    go_enrichment(., 
                  org = sample.name[[i]], 
                  function.dir = path_to_repo,
                  bg = 'all') %>% # enrichment against all ophio_cflo genes in the genome
    
    # # pull gene names for a given GO term
    # separate_rows(., gene_name, sep = ", ") %>%
    # filter(GO == "GO:0009405") %>% # pathogenesis
    # # filter(GO == "GO:0090729") %>% # toxin activity
    # # filter(GO == "GO:0044419") %>% # interspecies interaction between organisms
    # # filter(GO == "GO:0020037") %>% # heme binding
    # pull()
    
    go_enrichment_plot(clean = "no",
                       function.dir = path_to_repo) %>% 
    print()
  
}

## running GO enrichment for NOT-EXPRESSED genes in beau
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 10364"
## [1] "Number of genes in the test set: 759"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2456"
## [1] "Number of GO terms (at least 5genes) in background geneset: 925"
## [1] "Number of GO terms (at least 5genes) in test set: 140"
## [1] "Testing for enrichment..."

## running GO enrichment for NOT-EXPRESSED genes in ophio_cflo
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 7455"
## [1] "Number of genes in the test set: 190"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2400"
## [1] "Number of GO terms (at least 5genes) in background geneset: 844"
## [1] "Number of GO terms (at least 5genes) in test set: 50"
## [1] "Testing for enrichment..."

## running GO enrichment for NOT-EXPRESSED genes in ophio_kim
## [1] "Loading annotation file for Ophiocordyceps kimflemingae"
## [1] "Done."
## [1] "Number of genes in background geneset: 8441"
## [1] "Number of genes in the test set: 111"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 1099"
## [1] "Number of GO terms (at least 5genes) in background geneset: 279"
## [1] "Number of GO terms (at least 5genes) in test set: 4"
## [1] "Testing for enrichment..."

# B: genes that are expressed (FPKM > 1 for at least one time point)
expressed <- list()
for (i in 1:length(sample.name)) {
  expressed[[i]] <- 
    tbl(data.db, paste0(sample.name[[i]],"_expressed_genes")) %>% 
    filter(expressed=="yes") %>% 
    collect() %>% 
    pull(gene_name) 
  
  # How many genes are expressed?
  writeLines(paste("n(EXPRESSED) in", sample.name[[i]], ":", length(expressed[[i]])))
}

## n(EXPRESSED) in beau : 9006
## n(EXPRESSED) in ophio_cflo : 6998
## n(EXPRESSED) in ophio_kim : 8150

2. Daily rhythms in gene expression

## Load all the rhythmic genesets 
## Note, ordered according to their p-value; highly rhythmic at the top.
#
# Choose period
period = '24'

## 
rhy <- list()
for (i in 1:2) {
  rhy[[i]] <-
    tbl(ejtk.db, paste0(sample.name[[i]],"_zscores_",period,'h')) %>%
    filter(GammaP < gamma.pval) %>%
    select(ID, GammaP) %>% collect() %>% arrange(GammaP) %>%
    select(ID) %>% pull()
  
  # How many genes are rythmic?
  writeLines(paste0("n(rhythmic-",period, "h) in ", sample.name[[i]], " : ", length(rhy[[i]])))
}

## n(rhythmic-24h) in beau : 1872
## n(rhythmic-24h) in ophio_cflo : 2285

Hierarchical clustering of rhy24

• perform hierarchical clustering of 24h-rhythmic genes into four clusters;
• plot time-course heatmaps for the clustered 24h-rhythmic geneset
• Identify the day-peaking and night-peaking clusters visually.

# Choose period
period = '24'

## 
rhy <- list()
for (i in 1:2) {
  rhy[[i]] <-
    tbl(ejtk.db, paste0(sample.name[[i]],"_zscores_",period,'h')) %>%
    filter(GammaP < gamma.pval) %>%
    select(ID, GammaP) %>% collect() %>% arrange(GammaP) %>%
    select(ID) %>% pull()
  
  # How many genes are rythmic?
  writeLines(paste0("n(rhythmic-",period, "h) in ", sample.name[[i]], " : ", length(rhy[[i]])))
}

## n(rhythmic-24h) in beau : 1872
## n(rhythmic-24h) in ophio_cflo : 2285

## initialise lists to hold input and output of the hierarchical clustering
zscore.dat <- list() # zscore data (input)
my_gene_col <- list() # cluster identity for each rhythmic gene (output)
rhy.heat <- list() # pheatmap that can be saved/plotted (output)

# specify number of clusters
n_clusters <- 4

## run clustering and plot
for (i in 1:2) {
  ## load zscore dataset
  zscore.dat[[i]] <- data.db %>% tbl(., paste0(sample.name[[i]],"_zscores")) %>% collect()
  
  # Filter the zscores to keep only rhythmic genes
  zscore.rhy <-
    zscore.dat[[i]] %>% 
    filter(gene_name %in% rhy[[i]]) %>% 
    as.data.frame()
  
  # Set genes as rownames and convert it into a matrix
  rownames(zscore.rhy) = zscore.rhy$gene_name
  zscore.rhy <- as.matrix(zscore.rhy[-1])
  
  
  # Hierarchical clustering of the genesets
  my_hclust_gene <- hclust(dist(zscore.rhy), method = "complete")
  
  
  # Make annotations for the heatmaps
  my_clusters <- cutree(tree = as.dendrogram(my_hclust_gene), k = n_clusters) # k=  clusters
  my_gene_col[[i]] <- data.frame(cluster = my_clusters)
  
  
  # I’ll add some column annotations and create the heatmap.
  # Annotations for:
  # 1. Is the sample collected during the light or dark phase? 
  my_sample_col <- data.frame(phase = rep(c("light", "dark", "light"), c(5,6,1)))
  row.names(my_sample_col) <- colnames(zscore.rhy)
  
  # Manual color palette
  my_colour = list(
    phase = c(light = "#F2E205", dark = "#010440"),
    cluster = viridis::cividis(100)[c(10,90,60,30)]) #### NEED TO CHANGE #### account for n_clusters
  
  # Color scale
  my.breaks = seq(min(zscore.rhy), max(zscore.rhy), by=0.1)
  # my.breaks = seq(min(zscore.rhy), max(zscore.rhy), by=0.06)
  
  # Let's plot!
  rhy.heat[[i]] <-
    pheatmap(zscore.rhy, show_rownames = F, show_colnames = F,
             annotation_row = my_gene_col[[i]], 
             annotation_col = my_sample_col,
             cutree_rows = n_clusters, # OG was 4
             cutree_cols = 2,
             annotation_colors = my_colour,
             border_color=FALSE,
             cluster_cols = F,
             breaks = my.breaks,
             ## color scheme borrowed from: 
             color = inferno(length(my.breaks) - 1),
             # treeheight_row = 0, 
             # treeheight_col = 0,
             # remove the color scale or not
             main = paste0(sample.name[[i]], " 24h-rhythmic \n (n=", nrow(zscore.rhy), " genes)"),
             ## annotation legend
             annotation_legend = T,
             ## Color scale
             legend = T)
  
}

Phase plots

rhy.24.sig <- list()
phase.ejtk <- list()

# Obtain the phases of 24h-rhythmic genes beau v. ophio_cflo
for (i in 1:2) {

rhy.24.sig[[i]] <- 
  tbl(ejtk.db, paste0(sample.name[i],"_zscores_24h")) %>% 
  filter(GammaP < gamma.pval) %>% 
  collect()

# Get the phases of the best matched waveforms
phase.ejtk[[i]] <- circular::circular(rhy.24.sig[[i]]$Phase, units="hours", template="clock24")
# # Get the time-of-day of expression peak
# phase.ejtk[[i]] <- circular::circular(rhy.24.sig[[i]]$MaxLoc, units="hours", template="clock24")
# # Get the time-of-day of expression trough
# phase.ejtk[[i]] <- circular::circular(rhy.24.sig[[i]]$MinLoc, units="hours", template="clock24")

}

# save all the circular phases in a list
l.phases <- phase.ejtk
# let's name the list elements for later use and reference
names(l.phases) <- sample.name[1:2]

writeLines("Performing Watson test to check if the average peak of 24h-rhythms in Beau and Ophio-cflo differs significantly")

## Performing Watson test to check if the average peak of 24h-rhythms in Beau and Ophio-cflo differs significantly

# For all rhy genes
beau.ophio <- watson.two.test(l.phases[[1]],l.phases[[2]], alpha = FDR/100)
writeLines("Beau v. Ophio-cflo")

## Beau v. Ophio-cflo

beau.ophio %>% print()

## 
##       Watson's Two-Sample Test of Homogeneity 
## 
## Test Statistic: 5.9634 
## Level 0.05 Critical Value: 0.187 
## Reject Null Hypothesis

## Plot the phase distributions

# Initialize a list for saving the ggplots
g <- list()

means <- as.numeric(lapply(phase.ejtk, mean))
means <- circular(means, units="hours", template="clock24")


for(i in 1:length(l.phases)) {
  
  # define phase levels
  ordered_phases <- c("2","4","6","8","10","12",
                      "14","16","18","20","22","24")
    
  df.test <- l.phases[[i]] %>%
    as.data.frame() %>% 
    mutate(phase = x) %>%
    mutate(phase = replace(phase, x=="0", "24")) %>% 
    select(-x) %>% 
    group_by(phase = factor(phase, levels = ordered_phases)) %>%
    summarise(n_genes = n())

  m <- as.numeric(means[i])
    
  g[[i]] <-
    ggplot(df.test, aes(x=factor(phase), y=n_genes)) + 
    geom_bar(stat='identity', fill=col.scheme[[i]]) +
    
    xlab(c(names(l.phases)[i])) +
    
    scale_y_continuous(breaks = c(0,200,400,600)) +
    
    coord_polar() +
    theme_Publication() +
    theme(text = element_text(size = 15, colour = 'black'),
              # axis.title.x=element_blank(),
              # axis.text.x=element_blank(),
              legend.position = "none")
    #ggtitle(paste0("Dataset: ", names(l.phases)[i])) 
  
}

ggpubr::ggarrange(plotlist=g,
                  nrow = 2, ncol = 1,
                  widths = c(1,2), labels = NA)

Clusters of rhythmic genes

• Which processes are identified clusters of rhythmic genes involved in?
• How do they fluctuate throughout the day?

Beau - 24h

for (i in 1:n_clusters){
  
  writeLines(paste0("Species: ", sample.name[[1]], "\n", "24h-rhythmic genes, Cluster: ", i))
  
  # Summary
  genes <- my_gene_col[[1]] %>% rownames_to_column("g") %>% filter(cluster==as.character(i)) %>% pull(g)
  writeLines(paste0("n(genes) = ", length(genes),"\n"))
  
  # Enrichment
  overrepresented.terms <-
    genes %>% 
      go_enrichment(.,
                    function.dir = path_to_repo,
                    org = sample.name[[1]],
                    bg = expressed[[1]]) %>% 
      filter(adj_pVal < FDR/100) %>% 
      filter(over_under == "over")
  writeLines(paste0("\n", "n(overrepresented terms) = ", nrow(overrepresented.terms), "\n"))
  
  # Enriched terms word-cloud (borrowed from: https://towardsdatascience.com/create-a-word-cloud-with-r-bde3e7422e8a)
  if (nrow(overrepresented.terms)>0){
    # load libraries
    pacman::p_load(tm, wordcloud, RColorBrewer, wordcloud2)
    # get text as a character vector
    text <- overrepresented.terms %>% pull(GO_desc)
    # load your text data as a corpus
    docs <- Corpus(VectorSource(text)) # requires library "tm"
    # clean text (necessary?)
    docs <- docs %>%
              tm_map(removeNumbers) %>%
              tm_map(removePunctuation) %>%
              tm_map(stripWhitespace)
    docs <- tm_map(docs, content_transformer(tolower))
    docs <- tm_map(docs, removeWords, c("process", "molecular","cellular",
                                        "component", "compound", "part",
                                        "activity", "acid"
                                        ))
    # create document-term-matrix
    dtm <- TermDocumentMatrix(docs) 
    matrix <- as.matrix(dtm) 
    words <- sort(rowSums(matrix),decreasing=TRUE) 
    df <- data.frame(word = names(words),freq=words)
    # generate word-cloud
    wordcloud::wordcloud(words = df$word, freq = df$freq, min.freq = 2,
              max.words=200, random.order=FALSE, rot.per=0.35,
              scale=c(5,0.15),
              # colors=brewer.pal(8, "Dark2")
              colors=col.scheme[[1]]
              )
    
    ## save overrepresented GO terms for REVIGO analyses
      overrepresented.terms %>% select(GO, adj_pVal) %>% 
        readr::write_tsv(., paste0(path_to_repo,"/results/go_temp_files/",sample.name[[1]],"_Cluster_",i,".txt"))
      
  }
  # Stacked zplot
  genes %>% 
  stacked.zplot_tc6(cond = "beau") %>% 
    multi.plot(rows = 1, cols = 1)
  
}

## Species: beau
## 24h-rhythmic genes, Cluster: 1
## n(genes) = 767
## 
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 9006"
## [1] "Number of genes in the test set: 767"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2453"
## [1] "Number of GO terms (at least 5genes) in background geneset: 898"
## [1] "Number of GO terms (at least 5genes) in test set: 235"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 143

## Species: beau
## 24h-rhythmic genes, Cluster: 2
## n(genes) = 550
## 
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 9006"
## [1] "Number of genes in the test set: 550"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2453"
## [1] "Number of GO terms (at least 5genes) in background geneset: 898"
## [1] "Number of GO terms (at least 5genes) in test set: 179"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 2

## Species: beau
## 24h-rhythmic genes, Cluster: 3
## n(genes) = 337
## 
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 9006"
## [1] "Number of genes in the test set: 337"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2453"
## [1] "Number of GO terms (at least 5genes) in background geneset: 898"
## [1] "Number of GO terms (at least 5genes) in test set: 144"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 0

## Species: beau
## 24h-rhythmic genes, Cluster: 4
## n(genes) = 218
## 
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 9006"
## [1] "Number of genes in the test set: 218"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2453"
## [1] "Number of GO terms (at least 5genes) in background geneset: 898"
## [1] "Number of GO terms (at least 5genes) in test set: 122"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 36

Ophio-cflo - 24h

for (i in 1:n_clusters){
  
  writeLines(paste0("Species: ", sample.name[[2]], "\n", "24h-rhythmic genes, Cluster: ", i))
  
  # Summary
  genes <- my_gene_col[[2]] %>% rownames_to_column("g") %>% filter(cluster==as.character(i)) %>% pull(g)
  writeLines(paste0("n(genes) = ", length(genes),"\n"))
  
  # Enrichment
  overrepresented.terms <-
    genes %>% 
      go_enrichment(.,
                    function.dir = path_to_repo,
                    org = sample.name[[2]],
                    bg = expressed[[2]]) %>% 
      filter(adj_pVal < FDR/100) %>% 
      filter(over_under == "over")
  writeLines(paste0("\n", "n(overrepresented terms) = ", nrow(overrepresented.terms), "\n"))
  
  # Enriched terms word-cloud (borrowed from: https://towardsdatascience.com/create-a-word-cloud-with-r-bde3e7422e8a)
    if (nrow(overrepresented.terms) > 0) {
      # load libraries
      pacman::p_load(tm, wordcloud, RColorBrewer, wordcloud2)
      # get text as a character vector
      text <- overrepresented.terms %>% pull(GO_desc)
      # load your text data as a corpus
      docs <- Corpus(VectorSource(text)) # requires library "tm"
      # clean text (necessary?)
      docs <- docs %>%
                tm_map(removeNumbers) %>%
                tm_map(removePunctuation) %>%
                tm_map(stripWhitespace)
      docs <- tm_map(docs, content_transformer(tolower))
      docs <- tm_map(docs, removeWords, c("process", "molecular","cellular",
                                          "component", "compound", "part",
                                          "activity", "acid"
                                          ))
      # create document-term-matrix
      dtm <- TermDocumentMatrix(docs) 
      matrix <- as.matrix(dtm) 
      words <- sort(rowSums(matrix),decreasing=TRUE) 
      df <- data.frame(word = names(words),freq=words)
      # generate word-cloud
      wordcloud::wordcloud(words = df$word, freq = df$freq, min.freq = 2,
                max.words=200, random.order=FALSE, rot.per=0.35,
                scale=c(4,0.15),
                # colors=brewer.pal(8, "Dark2")
                colors=col.scheme[[2]]
                )
      
      ## save overrepresented GO terms for REVIGO analyses
      overrepresented.terms %>% select(GO, adj_pVal) %>% 
        # readr::write_tsv(., paste0("./results/go_temp_files/",sample.name[[2]],"_Cluster_",i,".txt"))
        write.table(., paste0(path_to_repo,"/results/go_temp_files/",sample.name[[2]],"_Cluster_",i,".csv"),
                    sep = ",",
                    row.names = F, col.names = F)
      
    }
  
  
  # Stacked zplot
  genes %>% 
  stacked.zplot_tc6(cond = "ophio_cflo") %>% 
    multi.plot(rows = 1, cols = 1)
  
}

## Species: ophio_cflo
## 24h-rhythmic genes, Cluster: 1
## n(genes) = 833
## 
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 6998"
## [1] "Number of genes in the test set: 833"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2359"
## [1] "Number of GO terms (at least 5genes) in background geneset: 831"
## [1] "Number of GO terms (at least 5genes) in test set: 221"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 83

## Species: ophio_cflo
## 24h-rhythmic genes, Cluster: 2
## n(genes) = 465
## 
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 6998"
## [1] "Number of genes in the test set: 465"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2359"
## [1] "Number of GO terms (at least 5genes) in background geneset: 831"
## [1] "Number of GO terms (at least 5genes) in test set: 158"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 31

## Species: ophio_cflo
## 24h-rhythmic genes, Cluster: 3
## n(genes) = 354
## 
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 6998"
## [1] "Number of genes in the test set: 354"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2359"
## [1] "Number of GO terms (at least 5genes) in background geneset: 831"
## [1] "Number of GO terms (at least 5genes) in test set: 158"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 30

## Species: ophio_cflo
## 24h-rhythmic genes, Cluster: 4
## n(genes) = 633
## 
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 6998"
## [1] "Number of genes in the test set: 633"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2359"
## [1] "Number of GO terms (at least 5genes) in background geneset: 831"
## [1] "Number of GO terms (at least 5genes) in test set: 242"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 123

Reduce redundant GO terms

Two options to do this:

Option #1: Get the list of overrepresented GO terms and their associated p-values and use REVIGO portal online to reduce the redundant terms

Option #2: Use the scripts provided by REVIGO to programmatically run REVIGO using bash/R. For more information see (here)[http://revigo.irb.hr/FAQ.aspx#q07] Status: tried running it via bash, and it didn’t work; NEED TO FIGURE IT OUT.

3. Beau v. Ophio - 24h-rhy

Next, we compare the homologous genes in both the fungi to understand if the rhythmic genes (and processes) in the two fungi are similar or not; also, is there any differences in the daily expression of these genes between the two fungal parasites?

Obtain homology data

# Read the source file
homology.file <- "ophio_beau_homology.csv"
homology.file <- 
  paste0(path_to_repo, "/results/proteinortho/", homology.file) %>% 
  read.csv(., stringsAsFactors = F, na.strings = c(" ","","NA"))

# Clean the source file to keep distinct gene-gene homologs
homology.dat <-
  homology.file %>% 
  # names() %>% 
  select(ophio_gene, beau_gene) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(beau_gene) %>% 
  filter(n()==1) %>% 
  select(beau_gene, ophio_gene)

writeLines(paste("Of the", length(expressed[[2]]), "genes expressed in Ophio-cflo,",
                 "and", length(expressed[[1]]), "genes expressed in Beau",
                 nrow(homology.dat), "genes show one-to-one orthology"))

## Of the 6998 genes expressed in Ophio-cflo, and 9006 genes expressed in Beau 5274 genes show one-to-one orthology

Summarize the results

for (i in 1:2){
  
  # exp.dat <- expressed[[i]]
  rhy.dat <- rhy[[i]]
  ortho.dat <- homology.dat %>% pull(i) 

  listInput <- list(rhy.dat, ortho.dat)
  names(listInput) <- c(paste0(sample.name[[i]], c("_rhy24","_ortho")))
  
  library(UpSetR)
  library(viridis)
  # caste.col <- c("#F23030","#1A80D9")
  upset(fromList(listInput), 
        number.angles = 0, point.size = 3, line.size = 1.5, 
        mainbar.y.label = "Number of overlapping genes", 
        sets.x.label = "Sig. rhy genes", 
        text.scale = c(1.5, # y-axis label ("# overlapping genes")
                       2, # y-axis tick labels ("1000, 2000,..")
                       1.5, # label for histogram ("sig. rhy genes")
                       1, # tick labels for histogram
                       1.5, # set names ("Cflo-brain_08h,..") 
                       1.5),
        sets = names(listInput),
        nintersects = 15,
        keep.order = T,
        sets.bar.color = viridis(1),
        # adding queries
        query.legend = "bottom"
        ) %>% 
    print()
}

Rhy24 genes w/ orthologs

rhy.homology.dat <- 
  homology.dat %>% 
  filter(beau_gene %in% rhy[[1]] | ophio_gene %in% rhy[[2]])

Hierarchical clustering

### Make the dataframe for plotting
zscore.rhy.homology.dat <-
  zscore.dat[[1]] %>% 
  filter(gene_name %in% rhy.homology.dat[[1]]) %>% 
  rename_at(vars(starts_with("ZT")), ~ (gsub("A", "B", .x, fixed = TRUE))) %>% # fix colnames for beau
  # add ophio homologs for the beau genes
  left_join(rhy.homology.dat, by=c("gene_name" = "beau_gene")) %>%
  # remove the beau names and keep the ophio names only
  select(-1) %>% 
  select(gene_name = ophio_gene, everything()) %>% 
  # join ophio-cflo data
  left_join(zscore.dat[[2]], by="gene_name") %>%
  # drop any genes without expression values (NA)
  na.omit() %>% 
  as.data.frame() %>% 
  # set genes as rownames
  column_to_rownames("gene_name")

# Set genes as rownames and convert it into a matrix
# rownames(zscore.rhy.homology.dat) = zscore.rhy.homology.dat$gene_name
zscore.rhy.homology.dat <- as.matrix(zscore.rhy.homology.dat)


# Hierarchical clustering of the genesets
my_hclust_gene <- hclust(dist(zscore.rhy.homology.dat), method = "complete")


# Make annotations for the heatmaps
n_clusters <- 4
my_clusters <- cutree(tree = as.dendrogram(my_hclust_gene), k = n_clusters) # k=  clusters
my_gene_col <- data.frame(cluster = my_clusters)


# I’ll add some column annotations and create the heatmap.
# Annotations for:
# 1. Is the sample collected during the light or dark phase? 
my_sample_col <- data.frame(phase = rep(rep(c("light", "dark", "light"),c(5,6,1)),2),
                            conds = rep(c("beau", "ophio_cflo"), each=12))
row.names(my_sample_col) <- colnames(zscore.rhy.homology.dat)

# Manual color palette
my_colour = list(
  phase = c(light = "#F2E205", dark = "#010440"),
  conds = c(beau = "#5A829F", ophio_cflo = "#AD212F"),
  cluster = viridis::cividis(100)[c(10,90,60,30)]) #### NEED TO CHANGE #### account for n_clusters

# Color scale
my.breaks = seq(min(zscore.rhy.homology.dat), max(zscore.rhy.homology.dat), by=0.1)
# my.breaks = seq(min(zscore.rhy), max(zscore.rhy), by=0.06)

# Let's plot!
pheatmap(zscore.rhy.homology.dat, show_rownames = F, show_colnames = F,
           annotation_row = my_gene_col, 
           annotation_col = my_sample_col,
           cutree_rows = n_clusters, # OG was 4
           cutree_cols = 4,
           annotation_colors = my_colour,
           border_color=FALSE,
           cluster_cols = F,
           breaks = my.breaks,
           ## color scheme borrowed from: 
           color = inferno(length(my.breaks) - 1),
           # treeheight_row = 0, 
           # treeheight_col = 0,
           # remove the color scale or not
           main = paste0("24h-rhythmic \n (n=",
                         nrow(zscore.rhy.homology.dat), " orthologous genes)"),
           ## annotation legend
           annotation_legend = T,
           ## Color scale
           legend = T)

Individual clusters

for (j in 1:2) {
  
  for (i in 1:n_clusters){
    
    writeLines(paste0("Species: ", sample.name[[j]], "\n", "24h-rhythmic genes, Cluster: ", i))
    
    # Summary
    genes <- my_gene_col %>% rownames_to_column("g") %>% filter(cluster==as.character(i)) %>% pull(g)
    writeLines(paste0("n(genes) = ", length(genes),"\n"))
    
    # define the background geneset for enrichment analysis
    bg.genes <- homology.dat %>% pull(ophio_gene) %>% unique()
    
    ## Transform gene names (ophio -> beau) and refine background geneset
    if (j == 1) {
      genes <- 
        homology.dat %>% 
        filter(ophio_gene %in% genes) %>% 
        pull(beau_gene)
      bg.genes <- homology.dat %>% pull(beau_gene) %>% unique()
    }
    
    # Enrichment
    overrepresented.terms <-
      genes %>% 
        go_enrichment(.,
                      function.dir = path_to_repo,
                      org = sample.name[[j]],
                      bg = expressed[[j]]) %>% 
        filter(adj_pVal < FDR/100) %>% 
        filter(over_under == "over")
    writeLines(paste0("\n", "n(overrepresented terms) = ", nrow(overrepresented.terms), "\n"))
    
    # Enriched terms word-cloud (borrowed from: https://towardsdatascience.com/create-a-word-cloud-with-r-bde3e7422e8a)
    if (nrow(overrepresented.terms)>0){
      # load libraries
      pacman::p_load(tm, wordcloud, RColorBrewer, wordcloud2)
      # get text as a character vector
      text <- overrepresented.terms %>% pull(GO_desc)
      # load your text data as a corpus
      docs <- Corpus(VectorSource(text)) # requires library "tm"
      # clean text (necessary?)
      docs <- docs %>%
                tm_map(removeNumbers) %>%
                tm_map(removePunctuation) %>%
                tm_map(stripWhitespace)
      docs <- tm_map(docs, content_transformer(tolower))
      docs <- tm_map(docs, removeWords, c("process", "molecular","cellular",
                                          "component", "compound", "part",
                                          "activity", "acid"
                                          ))
      # create document-term-matrix
      dtm <- TermDocumentMatrix(docs) 
      matrix <- as.matrix(dtm) 
      words <- sort(rowSums(matrix),decreasing=TRUE) 
      df <- data.frame(word = names(words),freq=words)
      # generate word-cloud
      wordcloud::wordcloud(words = df$word, freq = df$freq, min.freq = 2,
                max.words=200, random.order=FALSE, rot.per=0,
                scale=c(4,0.15),
                # colors=brewer.pal(8, "Dark2")
                colors=col.scheme[[j]]
                )
      
      # ## save overrepresented GO terms for REVIGO analyses
      #   overrepresented.terms %>% select(GO, adj_pVal) %>% 
      #     readr::write_tsv(., paste0(path_to_repo,"/results/go_temp_files/",
      #                                sample.name[[1]],"_Cluster_",i,".txt"))
        
    }
    # Stacked zplot
    genes %>% 
    stacked.zplot_tc6(cond = sample.name[[j]]) %>% 
      multi.plot(rows = 1, cols = 1)
    
  }
}

## Species: beau
## 24h-rhythmic genes, Cluster: 1
## n(genes) = 718
## 
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 9006"
## [1] "Number of genes in the test set: 718"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2453"
## [1] "Number of GO terms (at least 5genes) in background geneset: 898"
## [1] "Number of GO terms (at least 5genes) in test set: 267"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 154

## Species: beau
## 24h-rhythmic genes, Cluster: 2
## n(genes) = 1026
## 
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 9006"
## [1] "Number of genes in the test set: 1026"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2453"
## [1] "Number of GO terms (at least 5genes) in background geneset: 898"
## [1] "Number of GO terms (at least 5genes) in test set: 263"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 134

## Species: beau
## 24h-rhythmic genes, Cluster: 3
## n(genes) = 451
## 
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 9006"
## [1] "Number of genes in the test set: 451"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2453"
## [1] "Number of GO terms (at least 5genes) in background geneset: 898"
## [1] "Number of GO terms (at least 5genes) in test set: 164"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 16

## Species: beau
## 24h-rhythmic genes, Cluster: 4
## n(genes) = 324
## 
## [1] "Loading annotation file for Beauveria bassiana"
## [1] "Done."
## [1] "Number of genes in background geneset: 9006"
## [1] "Number of genes in the test set: 324"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2453"
## [1] "Number of GO terms (at least 5genes) in background geneset: 898"
## [1] "Number of GO terms (at least 5genes) in test set: 177"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 95

## Species: ophio_cflo
## 24h-rhythmic genes, Cluster: 1
## n(genes) = 718
## 
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 6998"
## [1] "Number of genes in the test set: 718"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2359"
## [1] "Number of GO terms (at least 5genes) in background geneset: 831"
## [1] "Number of GO terms (at least 5genes) in test set: 269"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 153

## Species: ophio_cflo
## 24h-rhythmic genes, Cluster: 2
## n(genes) = 1026
## 
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 6998"
## [1] "Number of genes in the test set: 1026"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2359"
## [1] "Number of GO terms (at least 5genes) in background geneset: 831"
## [1] "Number of GO terms (at least 5genes) in test set: 269"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 131

## Species: ophio_cflo
## 24h-rhythmic genes, Cluster: 3
## n(genes) = 451
## 
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 6998"
## [1] "Number of genes in the test set: 451"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2359"
## [1] "Number of GO terms (at least 5genes) in background geneset: 831"
## [1] "Number of GO terms (at least 5genes) in test set: 162"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 35

## Species: ophio_cflo
## 24h-rhythmic genes, Cluster: 4
## n(genes) = 324
## 
## [1] "Loading annotation file for Ophiocordyceps camponoti-floridani"
## [1] "Done."
## [1] "Number of genes in background geneset: 6998"
## [1] "Number of genes in the test set: 324"
## [1] "--------------------------------"
## [1] "Number of GO terms in background geneset: 2359"
## [1] "Number of GO terms (at least 5genes) in background geneset: 831"
## [1] "Number of GO terms (at least 5genes) in test set: 176"
## [1] "Testing for enrichment..."
## 
## n(overrepresented terms) = 77

Overlap b/n clusters and rhy24

## Visualize the overlap

cluster.dat <- list()
for (i in 1:n_clusters) {
    cluster.dat[[i]] <- my_gene_col %>% 
rownames_to_column("g") %>% filter(cluster==as.character(i)) %>% pull(g)
}
names(cluster.dat) <- paste0("Cluster_",1:4)

for (j in 1:2) {
    
    rhy.dat <- rhy[[j]]
    cluster.dat.dummy <- cluster.dat
    
    if (j == 1) {
      for (i in 1:n_clusters) {
        cluster.dat.dummy[[i]] <- 
            homology.dat %>% 
            filter(ophio_gene %in% cluster.dat.dummy[[i]]) %>% 
            pull(beau_gene)
      }
    }
    
    listInput <- list(rhy.dat, 
                      cluster.dat.dummy[[1]], cluster.dat.dummy[[2]], 
                      cluster.dat.dummy[[3]], cluster.dat.dummy[[4]])
    names(listInput) <- c(paste0(sample.name[[j]], c("_rhy24")), paste0("cluster_",1:4))
    
    library(UpSetR)
    library(viridis)
    # caste.col <- c("#F23030","#1A80D9")
    upset(fromList(listInput), 
          number.angles = 0, point.size = 3, line.size = 1.5, 
          mainbar.y.label = "Number of overlapping genes", 
          sets.x.label = "Sig. rhy genes", 
          text.scale = c(1.5, # y-axis label ("# overlapping genes")
                         2, # y-axis tick labels ("1000, 2000,..")
                         1.5, # label for histogram ("sig. rhy genes")
                         1, # tick labels for histogram
                         1.5, # set names ("Cflo-brain_08h,..") 
                         1.5),
          sets = names(listInput),
          nintersects = 15,
          keep.order = T,
          sets.bar.color = viridis(1),
          # adding queries
          query.legend = "bottom"
          ) %>% 
      print()
  
}

It seems that majority of the genes in Cluster 3 and 4 are sig. rhythmic in Ophio but not in Beau. We will perform the pairwise Fisher’s exact test to find out. Let’s dig in!

NOTE: We need to think about the best way to perform the Fisher’s exact test. For starters, I am transforming all the gene names to Ophio

# LIST ONE - Cluster identity
list1 <- cluster.dat
names(list1) <- names(cluster.dat)

## LIST TWO - ophio rhythmic genes
beau.ortho.rhy <- homology.dat %>% filter(beau_gene %in% rhy[[1]]) %>% pull(ophio_gene) %>% unique() 
ocflo.ortho.rhy <- homology.dat %>% filter(ophio_gene %in% rhy[[2]]) %>% pull(ophio_gene) %>% unique()
list2 <- list(beau.ortho.rhy, ocflo.ortho.rhy)
names(list2) <- paste0(sample.name[1:2], c("_24h"))

## CHECK FOR OVERLAP
library(GeneOverlap)

# define the background geneset 
# in our case, it would be the number of orthologous genes between beau and Ophio_cflo
nGenes = homology.dat %>% pull(ophio_gene) %>% unique() %>% length()

## make a GOM object
gom.1v2 <- newGOM(list1, list2, genome.size = nGenes)
png(paste0(path_to_repo, "/results/figures/BD/ocflo_beau_orthologs_rhy_overlap.png"),
    width = 15, height = 15, units = "cm", res = 300)
drawHeatmap(gom.1v2,
              adj.p=T,
              cutoff=0.01,
              what="odds.ratio",
              # what="Jaccard",
              log.scale = T,
              note.col = "grey60")
trash <- dev.off()

Orthologous rhy24 genes

As we predicted, Cluster 3 and 4 genes show a stronger overlap with 24h-rhythmic genes in O. cflo as comapred to Beauveria (as can be seen from both log2-odds ratio and the associated q-value). The signal is strongest for Cluster 4, so let’s see which genes are in this cluster.

Cluster-4 genes

## Get the annotation data
ocflo.annots <- read.csv(paste0(path_to_repo, "/data/ophio_cflo_TC6_data.csv"), stringsAsFactors = F)
ocflo.annots %<>% 
  as.tibble() %>%
  filter(expressed=="yes") %>%
  select(gene_name = gene_ID_ncbi, gene_ID_robin, blast_annot, GammaP_24h, GOs:ophio_kim_homolog)

## Subset the gene_names
ocflo.annots %>% 
  filter(gene_name %in% cluster.dat[[4]]) %>%
  
  # # plot the q-values
  # ggplot() + 
  # geom_density(aes(x=GammaP_24h)) +
  # labs(x="rhythmicity, 24h (q-value)") +
  # scale_x_continuous(breaks = c(0,0.05, 0.1, 0.5, 1)) +
  # theme_Publication()
  
  filter(!is.na(GOs)) %>% 
  view()


## Check these genes for other annotations (signalP, SSP, TMHMM)
# LIST ONE - Cluster identity
list1 <- cluster.dat
names(list1) <- names(cluster.dat)

## LIST TWO - ophio rhythmic genes
signalP <- ocflo.annots %>% filter(signalP == "yes") %>% pull(gene_name)
SSP <- ocflo.annots %>% filter(SSP == "yes") %>% pull(gene_name)
TMHMM <- ocflo.annots %>% filter(TMHMM == "yes") %>% pull(gene_name)
list2 <- list(signalP, SSP, TMHMM)
names(list2) <- paste0(sample.name[[2]], "-", c("signalP", "SSP", "TMHMM"))

## CHECK FOR OVERLAP
library(GeneOverlap)

# define the background geneset 
# in our case, it would be the number of orthologous genes between beau and Ophio_cflo
nGenes = ocflo.annots %>% nrow() 

## make a GOM object
gom <- newGOM(list1, list2, genome.size = nGenes)
png(paste0(path_to_repo, "/results/figures/BD/ocflo_beau_orthologs_annots_overlap.png"),
    width = 15, height = 15, units = "cm", res = 300)
drawHeatmap(gom,
              adj.p=T,
              cutoff=0.01,
              what="odds.ratio",
              # what="Jaccard",
              log.scale = T,
              note.col = "grey60")
trash <- dev.off()

Overlap of orthologous rhy24 gene clusters with additional annotations

4. Ophio DEGs during infection

Do prep for analyses

Prepare the functions, libraries required

# Let's load functions for running limorhyde
source(system.file('extdata', 'vignette_functions.R', package = 'limorhyde'))
# Let's load the libraries required for running Limorhyde
# library('annotate')
library('data.table')
library('foreach')
# library('GEOquery')
library('ggplot2')
library('knitr')
library('limma')
library('limorhyde')

conflict_prefer("union", "dplyr")

Format metadata

Create dataframe with metadata information for the different samples collected

sampleName <- c("ophio_cflo","ophio_ophio-infected")
short.name <- c("AC","AI") # AC = arb2-control, AI = arb2-infection
time.points <- c(2,4,6,8,10,12,14,16,18,20,22,24)
light.dark <- c(rep("light",times=5), rep("dark",times=6), rep("light", times=1))

meta <- data.frame(title = paste0(rep(sampleName, each=12),"_ZT",time.points),
                       sample = paste0(rep(time.points, times=2),rep(short.name, each=12)),
                       genotype = rep(sampleName, each=12),
                       time = rep(time.points, times=2),
                       cond = rep(sampleName, each=12),
                       LD = rep(light.dark, times=2),
                       stringsAsFactors = F)

meta %>% glimpse()

## Observations: 24
## Variables: 6
## $ title    <chr> "ophio_cflo_ZT2", "ophio_cflo_ZT4", "ophio_cflo_ZT6", "ophio…
## $ sample   <chr> "2AC", "4AC", "6AC", "8AC", "10AC", "12AC", "14AC", "16AC", …
## $ genotype <chr> "ophio_cflo", "ophio_cflo", "ophio_cflo", "ophio_cflo", "oph…
## $ time     <dbl> 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 2, 4, 6, 8, 10, …
## $ cond     <chr> "ophio_cflo", "ophio_cflo", "ophio_cflo", "ophio_cflo", "oph…
## $ LD       <chr> "light", "light", "light", "light", "light", "dark", "dark",…

Now, format the metadata.

### 1.1.1 Format the meta-data ----------------
# load the meta-data
sm <- meta
# Let's format the columns in the right data-type
sm$time <- as.numeric(sm$time)
# sm$batch <- as.factor(sm$batch)
sm$LD <- as.factor(sm$LD)
# sm$location <- as.factor(sm$location)

# Let's get a glimpse of the metadata
sm %>% as_tibble() %>% head()

## # A tibble: 6 x 6
##   title           sample genotype    time cond       LD   
##   <chr>           <chr>  <chr>      <dbl> <chr>      <fct>
## 1 ophio_cflo_ZT2  2AC    ophio_cflo     2 ophio_cflo light
## 2 ophio_cflo_ZT4  4AC    ophio_cflo     4 ophio_cflo light
## 3 ophio_cflo_ZT6  6AC    ophio_cflo     6 ophio_cflo light
## 4 ophio_cflo_ZT8  8AC    ophio_cflo     8 ophio_cflo light
## 5 ophio_cflo_ZT10 10AC   ophio_cflo    10 ophio_cflo light
## 6 ophio_cflo_ZT12 12AC   ophio_cflo    12 ophio_cflo dark

# Next we use limorhyde to calculate time_cos and time_sin, which are based on the first
#harmonic of a Fourier decomposition of the time column, and append them to the sm data frame.
sm = cbind(sm, limorhyde(sm$time, 'time_'))
# convert the dataframe into a data.table
sm <- data.table(sm)
# check that it worked
sm[1:5, ]

##              title sample   genotype time       cond    LD      time_cos
## 1:  ophio_cflo_ZT2    2AC ophio_cflo    2 ophio_cflo light  8.660254e-01
## 2:  ophio_cflo_ZT4    4AC ophio_cflo    4 ophio_cflo light  5.000000e-01
## 3:  ophio_cflo_ZT6    6AC ophio_cflo    6 ophio_cflo light  6.123234e-17
## 4:  ophio_cflo_ZT8    8AC ophio_cflo    8 ophio_cflo light -5.000000e-01
## 5: ophio_cflo_ZT10   10AC ophio_cflo   10 ophio_cflo light -8.660254e-01
##     time_sin
## 1: 0.5000000
## 2: 0.8660254
## 3: 1.0000000
## 4: 0.8660254
## 5: 0.5000000

load data

## DATASET 1
## Load the control O.cflo data (from TC6)
ocflo.control.dat <-
  data.db %>% 
  tbl(., paste0(sampleName[[1]], "_fpkm")) %>% 
  select(gene_name, everything()) %>%
  collect()

## DATASET 2
## Load the O.cflo infection data from the mixed transcriptomics study (from TC7)
inf.db <- dbConnect(RSQLite::SQLite(),
                   paste0(path_to_repo,"/../Das_et_al_2022b/data/databases/TC7_data.db"))
# src_dbi(inf.db)
# extract the (gene-expr X time-point) data
ocflo.inf.dat <-
  inf.db %>%
  tbl(., paste0(sampleName[[2]], "_fpkm")) %>%
  select(gene_name, everything()) %>%
  collect()

filter data

The goal is to use only the genes that show expression (>1 FPKM) for at least half of the timepoints during the 24h day (i.e., 6 of the 12 timepoints).

## DATASET 1
n.exp.1 <- apply(ocflo.control.dat[-1], 1, function(x) sum(x>=1))
ocflo.control.dat <- ocflo.control.dat[which(n.exp.1>=6),]
colnames(ocflo.control.dat)[-1] <- paste0("ZT", meta[meta$cond==sampleName[[1]],] %>% pull(sample))

## DATASET 2
n.exp.2 <- apply(ocflo.inf.dat[-1], 1, function(x) sum(x>=1))
ocflo.inf.dat <- ocflo.inf.dat[which(n.exp.2>=6),]
colnames(ocflo.inf.dat)[-1] <- paste0("ZT", meta[meta$cond==sampleName[[2]],] %>% pull(sample))

## Use the genes that are expressed in both conditions
emat <-
  ocflo.control.dat %>% 
  filter(gene_name %in% ocflo.inf.dat$gene_name) %>% 
  left_join(ocflo.inf.dat, by="gene_name") %>% 
  as.data.frame()

Run analyses

Next, let’s perform the DEG analyses for the ophio-cflo (halfway through infection v. controls)

### Convert to a matrix

# save gene names as row names
rownames(emat) <- emat[,1]
emat <- emat[,-1]
# Need to make the emat into a matrix.
emat <- data.matrix(emat)
# log2 transform the data
emat <- log2(emat + 1)


### Set thresholds
# Set threshold for q-value and log2FC
q.threshold <- 0.05 # currently, using 5% FDR (BH adjusted p-value)
log2.foldchange <- 1 # thus, any gene with a 2^(log2.foldchange) fold change in it's expression

### Format the metadata, if necessary
# Filter the metadata according to your comparison
sm.sub <- sm %>% filter(cond %in% c(sampleName))
# Define the cond column as a factor
sm.sub$cond <- as.factor(sm.sub$cond)

### Let's run the DEG analyses
# Use the subsetted emat to find DEGs
design.deg = model.matrix(~ cond + time_cos + time_sin, data = sm.sub)
#
fit = lmFit(emat, design.deg)
fit = eBayes(fit, trend = TRUE)
# Take a look at the coefficients table
# fit$coefficients %>% head()
#
deLimma.deg = data.table(topTable(fit, coef = 2, number = Inf), keep.rownames = TRUE)
setnames(deLimma.deg, 'rn', 'gene_name')
deLimma.deg[, adj.P.Val := p.adjust(P.Value, method = 'BH')]
setorderv(deLimma.deg, 'adj.P.Val')

### Annotate the results
# Annotate the results to indicate the significant genes
all.DEGs <-
  deLimma.deg %>% 
  arrange(desc(abs(logFC)), adj.P.Val) %>% 
  mutate(sig = as.factor(ifelse(adj.P.Val < q.threshold & abs(logFC) >= log2.foldchange, "yes", "no"))) %>% 
  mutate(inf_v_control = as.factor(ifelse(sig=="yes", ifelse( logFC > 0, "up", "down" ), "NA"))) %>% 
  mutate(inf = sampleName[[2]])

### Summarize the results
writeLines(paste0("\nControl-", sampleName[[1]], " v. ", sampleName[[2]], "\n--Results of DEG analysis--"))

## 
## Control-ophio_cflo v. ophio_ophio-infected
## --Results of DEG analysis--

## How many DEGs - 5% FDR and ≥ 1 fold change in gene expression
all.DEGs %>% 
  # filter(adj.P.Val < q.threshold) %>%
  # filter(abs(logFC) >= 2) %>% # change the criteria here for top DEG or all DEG (logFC≥1)
  filter(sig == "yes") %>% 
  
  # pull(gene_name) %>% 
  
  group_by(inf_v_control) %>% 
  summarise(n_genes = n()) %>% 
  as.data.frame() %>% 
    ## n = 81 up- and 141 down-regulated genes in Cflo heads during Ophio-infection 
    ## (at 5% FDR; log2-fold-change ≥ 1) 
  print()

##   inf_v_control n_genes
## 1          down     395
## 2            up     318

### Subset to keep only sig. DEGs
sig.DEGs <- all.DEGs %>% filter(sig=="yes")

Visualize results

# Volcano plot 
library(viridis)

ggplot(all.DEGs) +
  # geom_hline(yintercept = -log10(0.05), col="red", alpha=0.6) +
  # geom_vline(xintercept = c(-2,2), col="grey60", alpha=0.75) +
  geom_point(aes(x = logFC, y = -log10(adj.P.Val), color=sig), size = 1.5, alpha = 0.5) +
  labs(x = expression(log[2]*' fold-change (inf_v_control)'), 
       y = expression(-log[10]*' '*q[DE]),
       title = "O.cflo (infection v. control)",
       color = "significant") +
  # scale_x_continuous(limits = c(-5,3),
  #                    breaks = c(-5,-4,-3,-2,-1,0,1,2,3),
  #                    labels = c("-5","","-3","","-1","","1","","3")) +
  # xlim(c(-50,50)) +
  theme_Publication() +
  scale_color_viridis(discrete = T, direction = -1, option = "viridis")

Load manip data

## Load the ophio DEG (at manipulation) data from Will et al. 2020
will2020_data <- read.csv(paste0(path_to_repo,"/data/input/ophio_cflo/complete_annotations/FullBlast_EC05_RNAseq_orignal_copy_26Aug19.csv"), stringsAsFactors = F)
will2020_data %<>% 
  as_tibble() %>% 
  filter(sample_1=="Alive" & sample_2=="Fungus") %>%
  select(arb2_gene, logFC = log2.fold_change., q_value, significant) %>% 
  mutate(logFC=as.numeric(logFC), q_value=as.numeric(q_value)) %>%
  filter(significant=="yes") %>% 
  mutate(up_down = ifelse(logFC > 0, "down", "up")) %>%
  mutate(logFC = -1*logFC) %>% 
  na.omit()

### Change ophio gene names to ncbi IDs
will2020_data %<>% 
  left_join(ocflo.annots[1:3], by=c("arb2_gene"="gene_ID_robin")) %>%
  select(-1) %>% 
  select(gene_name, blast_annot, everything())

5. Overlap between DEGs during infection v. manipulation

### Subset the up/down-regulated genes
### At halfway-through disease progression
inf.up <- sig.DEGs %>% filter(inf_v_control=="up") %>% pull(gene_name)
inf.down <- sig.DEGs %>% filter(inf_v_control=="down") %>% pull(gene_name)
### At active manipulation
manip.up <- will2020_data %>% filter(up_down=="up") %>% pull(gene_name)
manip.down <- will2020_data %>% filter(up_down=="down") %>% pull(gene_name)


### Visualize the results
listInput <- list(inf.up, inf.down, manip.up, manip.down)
names(listInput) <- c(paste0("inf_",c("up","down")), paste0("manip_", c("up","down")))
    
library(UpSetR)
library(viridis)
upset(fromList(listInput), 
  number.angles = 0, point.size = 3, line.size = 1.5, 
  mainbar.y.label = "Number of overlapping genes", 
  sets.x.label = "Sig. DE genes", 
  text.scale = c(1.5, # y-axis label ("# overlapping genes")
                 2, # y-axis tick labels ("1000, 2000,..")
                 1.5, # label for histogram ("sig. rhy genes")
                 1, # tick labels for histogram
                 1.5, # set names ("Cflo-brain_08h,..") 
                 1.5),
  sets = names(listInput),
  nintersects = 15,
  keep.order = T,
  sets.bar.color = viridis(1),
  # adding queries
  query.legend = "bottom"
  ) %>% 
print()

### Test significance of overlap
list1 <- list(inf.up, inf.down)
names(list1) <- paste0("inf_",c("up","down"))
list2 <- list(manip.up, manip.down)
names(list2) <- paste0("manip_", c("up","down"))
bg.genes <- all.DEGs %>% nrow()

overlap <- check_overlap(list1, list2, bg.genes)

## $rowInd
## [1] 2 1
## 
## $colInd
## [1] 1 2
## 
## $call
## heatmap.2(x = plot.mat, Rowv = NA, Colv = NA, dendrogram = "none", 
##     scale = "none", col = brewer.pal(ncolused, grid.col), colsep = col_sep, 
##     rowsep = row_sep, sepcolor = "white", sepwidth = c(0.002, 
##         0.002), cellnote = note.mat, notecex = 1.6, notecol = note.col, 
##     trace = "none", margins = margins_use, cexRow = row_cexrc, 
##     cexCol = col_cexrc, key = T, keysize = key_size, density.info = "none", 
##     main = main.txt, xlab = footnote)
## 
## $carpet
##            inf_down   inf_up
## manip_up    0.00000 2.075416
## manip_down  1.20292 0.000000
## 
## $rowDendrogram
## 'dendrogram' with 2 branches and 2 members total, at height 1.414214 
## 
## $colDendrogram
## 'dendrogram' with 2 branches and 2 members total, at height 1.414214 
## 
## $breaks
##  [1] 0.0000000 0.2306017 0.4612035 0.6918052 0.9224070 1.1530087 1.3836105
##  [8] 1.6142122 1.8448140 2.0754157
## 
## $col
## [1] "#F7FCF5" "#E5F5E0" "#C7E9C0" "#A1D99B" "#74C476" "#41AB5D" "#238B45"
## [8] "#006D2C" "#00441B"
## 
## $colorTable
##         low      high   color
## 1 0.0000000 0.2306017 #F7FCF5
## 2 0.2306017 0.4612035 #E5F5E0
## 3 0.4612035 0.6918052 #C7E9C0
## 4 0.6918052 0.9224070 #A1D99B
## 5 0.9224070 1.1530087 #74C476
## 6 1.1530087 1.3836105 #41AB5D
## 7 1.3836105 1.6142122 #238B45
## 8 1.6142122 1.8448140 #006D2C
## 9 1.8448140 2.0754157 #00441B
## 
## $layout
## $layout$lmat
##      [,1] [,2]
## [1,]    4    3
## [2,]    2    1
## 
## $layout$lhei
## [1] 1.485097 4.000000
## 
## $layout$lwid
## [1] 1.485097 4.000000